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ABSTRACT 

We use orbit based dynamical models to fit the 2nd and 4th moments of the line of sight velocity distributions of the 
Fornax, Sculptor, Carina and Sextans dwarf spheroidal galaxies. Our goal is to compare dark matter profile models of 
these four systems using Bayesian evidence. We consider NFW, Einasto and several cored profiles for their dark halos 
and present the probability distribution functions of the model parameters. When considering each system separately, 
we find there is no preference for one of these specific parametric density profiles. However, the combined evidence shows 
that it is unlikely that all galaxies are embedded in the same type of cored profiles of the form poM oc 1/(1 + r^)°'^, 
where /? = 3, 4. For each galaxy, we also obtain an almost model independent, and therefore accurate, measurement 
of the logarithmic slope of the dark matter density distribution at a radius ~ r_3, i.e. where the logarithmic slope of 
the stellar density profile is —3. This slope ranges from ^ —1.4 for Fornax to ^ —1.1 for Sextans, both at '--^ 1 kpc. 
All our best fit models essentially have the same mass distribution over a large range in radius (from just below r_3 to 
the last measured data point). This remarkable finding likely implies much stronger constraints on the characteristics 
that subhalos extracted from cosmological simulations should have in order to host the dSph galaxies around the Milky 
Way. 
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1. Introduction 
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'According to galaxy formation theories dwarf spheroidal 
galaxies are believed to inhabit massive dark matter halos. 
Because of their large mass to light ratio these galaxies are 
ideal to test fundamental predictions of the ACDM cosmo- 
logical paradigm, since it is generally considered relatively 
safe to neglect baryons in the construction of dynamical 
models. 

One of the strongest predictions from ACDM concerns 
the dark matter density profile. Early simulations of dark 
matter halos assembled in a cosmological context showed 
that such a profile is accurately desc ribed by a two-sloped 
form , now known as NFW profile (jNavarro at al.l Il996l . 

I'IQQTT ). More recently Eina sto profiles hav e been shown 
ibo pr ovide a better fit (e.g. iNavarro at al.1 2010; Springel 
et al., 
fet al. 



200E )■ in particular for satellite galaxies ( Vera-Ciro 



2013|). These predictions are made using dark mat- 



ter only simulations and therefore neglect (by construc- 
tion) the baryonic component. And although baryons are 
sub-do minant in the total potential of the system (Walkeij 
l2012af ) , it has been suggested that they could play a role in 
the evolution of dwarf spheroidal galaxi es, for instance, in 
modifying the internal orbital structure (iBrvan et al.ll2012 ) 
and the overall density profile (Gov ernato et al.ll2012D . The 
complex evolution of baryons and its non-trivial interplay 
with the host halo are difficult to model and not ye t com- 
pletely understood (see lPontzen fc Governatoll2012f ). 

Another effect driving the internal dynamics of satel- 
lite halos is the tidal interactio n with the ma i n hos t. It 
can change the density profile (jHavashi et al.l 120031 ). the 



geom etrical shape of the mass distribution (jKuhlen et al.l 
I2007D, and also infiuen ce the kinematics of the embedded 
stars (jLokas et al.ll2010l ). Unfortunately these uncertainties 
imply that even when the observations of the local dwarf 
spheroidal galaxies are not consistent with being embed- 
ded in the halos predicted from pure dark matter N-body 
simulations, this does not necessarily refiect a fundamental 
problem of ACDM. 

Thanks to their relative proximity, information for indi- 
vidual stars in the dwarf galaxies satellites of the Milky Way 
are relatively easy to get. Sky positions are easily deter- 
mined from photometry, and radial velocity measurements 
are possible to estimate within an error of ~ 2 km/s. Some 
of the datasets compiled to the date include thousands of 
individual members with line-of-sight velocities (Battaglia 
et al. 2008 . 2011 . 2006 ; Helmi et al. 2006 : Walker et al] 
l2009al ) Proper motions of individual stars are currently still 
too difficult to measure. Despite the fact that only three of 
the total of six phase space coordinates are available from 
measurements, it is possible to create dynamical models of 
these systems that can be compared to these observables. 

Follow ing the method thoroughlv described in Bred- 
dels et al. (|2012| ) we set out to model Fornax. Sculptor. 
Carina and Sextans with orbit-based dynamical methods 
(Schwarzschild modeling) assuming they are embedded in 
spherical halos. As extensiv ely shown in the literature (e.g. 
ICretton et allllQQgl: I Jardel fc Gcbhardt 2012LRichstone & 
Tremaine'l984':'Rix et al.lll997tlValluri et al. 200^: van den 
Bosch et al. 2008; van der Marel et al.l ll998) this method 
allows to construct a non-parametric estimator of the dis- 
tribution function. Among many, this method has one ad- 
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Name 


A^Batt 


Mvalkcr 


-^*mcinbcr 


Fornax 


945(1) 


2633(^) 


2936 


Sculptor 


1073(2) 


1541(5) 


1685 


Carina 


811(3) 


1982(^) 


885 


Sextans 


792(4) 


947(5) 


541 



Table 1. Number of stars in the kinematic samples used in this 
paper. Sources: ^^^Battagli a et ali (|2006l ). ^^ iBattaglia et aU 
IpOOg) (^^Helm i et al. ( 20061): iKoch et all l(20M): Starkenburg 
et al. l|2010D ,^ {BattagUa et all (|201ir ). ^^ Malker eTal] (|2009al ) 



vantage over Jeans modeling, by not having to assume a 
particular velocity anisotropy profile, therefore being more 
general and thus less prone to biases associated to the as- 
sumptions. But even in this case there are other limitations 
in the modeling such as the mass-anisotropy degeneracy. In 
this work we use higher moments (fourth moment) of the 
line of sight velocity distribution to get a better handle on 
this degeneracy. 

To compare how different shapes for the dark matter 
profiles fit the data, we first need to establish a statis- 
tical framework. In this paper we do this in a Bayesian 
way using the evidence (Mackay 2003). This method pro- 
vides a natural way of comparing models in Bayesian in- 
ference and also makes it possible to combine the data of 
all the dwarf spheroidals to test for example, if all dwarf 
spheroidals could be embedded in a universal halo ( Mateo 
et al. 1993; Walker et al. 2009^). Furthermore, the shape 
may give us hints to how the dwarf galaxy formed and the 
anisotropy profile may be u sed to distinguish between evo- 
lutio narv scenarios (see e.g.lHelmi et al.l 20121: Kazantzidis 
et al. l2011l : lMavej|2010l ). 

This paper is organized as follows. We begin in ^ 
by presenting the data and all the ingredients needed to 
do the model comparison. In ^ we present our dynami- 
cal and statistical methods. We present the results of our 
Schwarzschild models for the four dSph in our sample in 
i j4.11 while the Bayesian model comparison is done in iJ4.2l 
We discuss the implications of our results in ij4.3l and con- 
clude in ^ 



2. Data 

In this section we present the data that is used for fitting 
our dynamical models. The radial velocit y measurements 
of the dwarf spheroidal galaxies come from|B attaglia et ahl 
( 20081. 120061) : iHelmi et all (|2006l 'l: IWalker etal. (2009ai ) and 
iBattaglia et al.1 ()2011l ). We plot radius versus heliocentric 
velocity in Fig. [T] for each galaxy separately. 

Figure [1] shows that each dSph suffers from foreground 
(Milky Way) contamination. To remove this contamina- 
tion and reliably identify member stars we have developed 
a simple analytic model for the positional and kinematic 
distribution of both foregroun d and the gal axy in ques- 
tion (along the lines of Batta glia et al.lj2008l . Breddels et 
al. in prep). For each particular dataselQ, we assume that 
the foreground has a constant surface density, and that the 
dSph follows a specific stellar density profile. We also as- 
sume that the velocity distribution at each radius may be 



modeled as sum of two Gaussians. The Gaussian describ- 
ing the foreground has the same shape at all radii, while 
that of the stars associated with the dwarf can have a vary- 
ing dispersion with radius. Their relative amplitude also 
changes as function of distance from the dwarf's centre. 
This model results in a determination of the relative con- 
tribution of member-to-non-member stars as a function of 
velocity and radial distance R. 

Based on this model we calculate the elliptical radius at 
which the ratio of dSph: foreground is 3:1 (without using any 
velocity information). We remove all stars outside of this 
radius from the dataset. A particular star included in more 
than one dataset is removed only when it it satisfies the 
condition for all sets, for instance a star outside i?e,cut, Batt, 
but inside i?e,cut, Walker wiU not be discarded. This simple 
clipping in elliptical radius cleans up part of the foreground 
contamination. For completeness, the radii for all datasets 
cleaned up in this way are presented in Table [21 as well as 
the fit to the foreground model. The number of stars and 
the sources are listed in Table [T] 

For the resulting dataset, we compute the second and 
fourth moment of the radial velocity as a function of cir- 
cular radius as followfl We first define radial bins such 
that each bin has at least 250 stars in the velocity range 
Wsys — 3CTt,,Wsys + 3cti, . If thc last bin has less that 150 ob- 
jects, the last two bins are merged. After this, we fit our 
parametric model for the galaxy plus foreground for each 
radial bin, to derive new velocity dispersions. Then for each 
bin we do a 3ct clipping on the velocity using the new ve- 
locity dispersion, and from this selection we calculate the 
second and fourth moments. The errors o n the moments 
are co mputed using Eqs. (17) and (19) in iBreddels et al.l 
(|2012f ). The second moment and the kurtosiqf] are shown 
in Fig. [2] for each galaxy, where the black dot corresponds 
to the mean, and the error bars indicate the la error bar. 
The blue region shows the confidence interval for the NFW 
model found in ii4.1l 

For the photometry we use analytic fits given by various 
literature sources as listed in Table [3] Although the stellar 
mass is sub-dominant in the potential, we do include its 
co ntribution in the dyna mic models and fix M/Ly = 1, as 
in IBreddels et al.l (|2012f ). 



3. Methods 

3.1. Dynamical models 

Our aim is to compare different models to establish what 
type of dark matter profile best matches the kinematical 
data of local dSph galaxies. Here we consider the following 
profiles to describe the dark matter halos of the dwarfs in 



^ For a given dSph there may be multiple datasets, and we treat 
each independently because their sampling might be diflerent. 



^ Elliptical radii are only used for the clipping, for the rest of 
the analysis we use the circular radius 

^ The kurtosis is defined as 72 = 114/ fi2, where fi4 is the fourth 
and /i2 is the second moment of the line of sight velocity distri- 
bution. 
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A. Helmi: Model comparison of dark matter profiles in local dSplis 
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Fig. 1. Radius versus line of sight velocity for Fornax, Sculptor Carina and Sextans. The horizontal lines show the borders of 
the bins, the vertical lines denote the mean systemic velocity of the galaxy together with the ±3(t region. 



Name 


-n^e.max,Batt 


J^e, max. walker 


Mmw 


(Tmw 


Mdwarf 


f^dwarf 




(kpc) 


(kpc) 


(km/s) 


(km/s) 


(km/s) 


(km/s) 


Fornax 


1.82 


2.21 


41.1 


38.9 


55.1 


12.1 


Sculptor 


1.37 


1.65 


17.9 


47.4 


110.6 


10.1 


Carina 


0.86 


0.96 


70.9 


62.5 


222.9 


6.6 


Sextans 


1.86 


1.65 


67.5 


74.5 


224.3 


7.9 



Table 2. Parameters of the foreground plus dwarf galaxy model used for determining membership, as well as for deriving the 
radial profiles for the second and fourth velocity moments for each dSph. 



Name 


distance 


profile 


scale radius 


Lv 




(kpc) 




(kpc) 


xIO^Lq 


Fornax 


138(1) 


Plummer^ 


0.79 


100(2) 


Sculptor 


79(3) 


Plummer-^ 


0.30 


10(3) 


Carina 


lOlW 


Exponential* 


0.16 


2.4(4) 


Sextans 


86(1) 


Exponential'! 


0.39 


4.37(4) 



Ta ble 3. Distan ces, type of photome t ric pr ofile used, scale radi u s and stell ar luminosity used for t h e dyn amic models. Sources: 
ri) lMated(|1998l ). ('2) iBattadia et al.l (|2006l ). fS) iBattadia et al.l (|2008l ). (A') llrwin fc Hatzidimitrioul (119951 ) 



our sample: 
p{r) 



Po 



x{l + x) 



2 ' 



p(r) = 



Po 



{1 + x^f^^ 
2 



p{r) = Po exp 



U-'-O) 



NFW (1) 

(cored) /37-profile (2) 

Einasto (3) 



where x = r/r^ and r^ is the scale radius. Each model has 
at least two unknown parameters r^ and po- As we did in 
[Breddels et al. (2012), we transform these two parameters 
to Tg and Mikpc (the mass within 1 kpc). As discussed in 
the Introduction, the NFW and Einasto models are known 



to fit the halos dark matter distributions extracted from 
cosmological N-body simulations. On the other hand, we 
explore the f3j models to test the possibility of a core in the 
dark halo. Note that, in comparison to the NFW profile, 
the Einasto model has one extra parameter (a'), but here 
we consider only t wo values for a' = 0.2, 0.4 to cover the 
range suggested by IVera-Ciro et aD (J2013l ). On the other 
hand, the f3j profiles have two extra parameters, but we 
limit ourselves here to two different outer slopes (/? — 3,4) 
and two different transition speeds between the inner and 
the outer slopes (7 — 1,2). Note that the /37 models have 
a true core only for 7 > 1, however in all cases the central 
logarithmic slope vanishes, dlogp/dlogr — 0. However, we 
loosely refer to these models as cored in what follows. Note 
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Fig. 2. Line of siglit velocity moments for Fornax, Sculptor Carina and Sextans. For eacli galaxy we show the velocity dispersion 
and the kurtosis. The blac k dots show the mean , and the error bars the la error. The blue regions show the confidence interval 
for the NFW fit, similar to lBreddels et al.l (|2012f ). 



that, with these choices, all of our profiles ultimately have 
just two free parameters. The list of models explored and 
their parameters are summarized in Table 21 

The orbit-based dynamical (Schwarzschild) models of 
each dwarf galaxy are obtaine d as fol lows (see Breddels 
et al. I2OI2I . for a more detailed description). For each of 
the dark halo profiles, with its own set of parameters, we 
integrate a large number of orbits in the respective grav- 
itational potential (including also the contribution of the 
stars). We then find a linear combination of these orbits 
that fits both the light and the kinematics. The orbital 
weights found in this way have a physical meaning and can 



be used to obtain the distribution function of the system. 
As data we have the line of sight velocity moments (sec- 
ond and fourth depicted in Fig. [5]), and the light profile 
(Table [5]) . The best fit models (which give us the values 
of the parameters for a specific dark matter halo profile) 
are those that minimize the x^ = Xkin + Xregi under the 
condition that the orbital weights are positive, and that 
the observed light distribution is fit to better than 1% at 
each radius. Here xlin = Efc(M2,fc - A<"fc''''0^/var(^2,fe) + 
Y.kif^4,k - M™fc''''0^/var(/X4^fe). The x^g is a regularization 
term to make sure that the solution for the orbit weights 
leads to a relatively smooth distribution function. Bred- 
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Name 


Fixed parameters 


Free parameters 


NFW 


- 


Mikpc, rs 


Corel 3 


/3 = 3,7 = 1 


Mikpc, rs 


corel4 


^ = 4,7 = 1 


Mikpc, Ts 


core23 


/3-3,7 = 2 


Mikpc, Ts 


core24 


/? = 4,7 = 2 


Mikpc, Ts 


einasto.2 


a' = 0.2 


Mikpc, Ts 


einasto.4 


a' = 0.4 


Mikpc, rs 



Table 4. Model names and their characteristic parameters of 
the variouR dark matter density profiles explored. 



dels et al. (120121 ) calibrated the amplitude of this term 
for Sculptor. To have the regularization term for the 
other dwarfs of the same relative strength, we note that 
Xreg Of 1/A^, where N is the number of members with radial 
velocities, since the x^jj^ term also scales as 1/iV. There- 
fore, normalizing its amplitude to that of Sculptor we get 

Xreg, dwarf " Xreg, Scl ^ -^Scl/^dwarf • 



3.2. Bayesian model comparison 

Background on B ayesian model comparison may be found 
in .Mackavl (|2003( ). For completeness we discuss it here 
briefly, but we assume the reader is familiar with the basics 
of Bayesian inference. 

Given the data V and assuming a model M^, the poste- 
rior for the parameters 9i of this model is: 



pi9,\V,M,) 



p{V%,Kh)p{e,\Mi 
p{V\M,) 



(4) 



The normalization constant p{'D\Mi), also called the ev- 
idence, is of little interest in parameter inference, but is 
useful in Bayesian model comparison. To assess the proba- 
bility of a particular model given the data, we find 



p{M,\V) = 



p{V\M,)p{Mi) 
p{V) ' 



(5) 



where we see the evidence is needed. In this case p{'D) is 
the uninteresting normalization constant, as it cancels out 
if we compare two models: 



p{M,\V) _ p{V\M,) p{M,) 
p{M,\V) p{V\M,)p{M,) 



B, 



'piM.y 



(6) 



where Bij is called the Bayes factor. If we take the priors 
on the different models to be equal (i.e. p{Mi) — p{Mj)), 
the ratio of the evidence (the Bayes factor Bij) gives the 
odds ratio of the two models given the data V. 

Using these results we can perform model compari- 
son between dark matter density profiles, i.e. Ai = 
{Mnfw, MEinasto, ■•■}, ^ud Calculate for instance the odds 
that a given galaxy is embedded in an NFW profile com- 
pared to an Einasto model, -BNFW.Einasto- 

Not only can we do model comparison on a single object, 
but we may also test if our objects share a particular model 
(e.g. they are all embedded in NFW halos). If our dataset 
V consists of the observations of two galaxies, i.e. V = 
Vi U 7)2 and assuming the datasets are uncorrelated and 
independent, we obtain: 



p(M,\V) _ p{Vi\M,)p{V2\Mi) piM,[ 



p{M, \V) p{V, |Mj)p(2?2 1 Af,) p{Mj 



— Bj^i^B, 



ij,l-Oij,2 



p{M,) 



(7) 

where each factor p{T>k\Mi) should be marginalized over its 
(own) characteristic parameters. From Eq. ([7]) we can see 
that the odds ratio of the models and Bayes factor from 
different measurements can be multiplied to give combined 
evidence for a particular model. 

Behind each p(Afi|2?) is a set of orbit based dynamical 
(Schwarzschild) models, obtained as described above. For 
each of the models we calculate the evidence. Later on we 
compare each model's evidence to that of an NFW profile, 
i.e. we compute the Bayes factor -Bi^NFW, where i can be e.g. 
Einasto. By definition i?NFW,NFW — 1, and again assuming 
equal priors on the different models, the Bayes factor equals 
the odds ratio of the models, such that for -Bi^NFW > 0, 
model i is favored over an NFW profile. 



4. Results 

4.1. Schwarzschild models 

As a result of our Schwarzschild modeling technique, we 
obtain a two dimensional probability density function (pdf ) 
of the two parameters, Mikpc and Tg, for each galaxy and 
for each dark matter halo profile. In Fig. [3] we plot the 
pdf for the cored models and the NFW and Einasto models 
separately for each galaxy. The colored dots correspond 
to the maximum likelihood for each of the corresponding 
models as indicated by the legend. The contours show the 
1 and 2cr equivalent confidence intervals (the 3cr contour is 
not shown for clarity). For both Fornax and Sculptor the 
parameters for all profiles are relatively well determined, 
while for Carina and Sextans this is less so. This can be 
attributed to the difference in sample size (and hence to 
the smaller number of members) in these systems, which 
has translated into fewer bins where the moments can be 
computed (see Fig. [5]). In general for all four galaxies the 
scale radius for the cored profiles is found to be smaller 
than that for the NFW/Einasto profiles. We come back to 
this point in section H75] 

Our model's masses at ri/2, the 3d radius en closing half 
of the stellar mass, are compatible with those of IWolf et al.l 
(|2010r ). However, o ur results for Fornax do not agree with 
those of iJardel fc G ebhardt (2012). These authors prefer a 
cored profile with a much larger scale r adius, and their en - 
closed mass is smaller in comparison to IWolf et aTl (|2010f ) . 
We note that this might be partly related to the fact that 
the amplitude of their line of sight velocity dispersion pro- 
file (see their Fig. 2) is slightly lower than what we have 
determined here. 

In Fig. [5] we overlay on the kinematic observables the 
predictions from the best fit Schwarzschild models. We note 
that all models provide very similar and virtually indistin- 
guishable fits, especially for the 2nd moment. Some slight 
differences are apparent in the kurtosis, but in all cases, the 
differences are smaller than the error bars on the moments. 

In general we find all anisotropy profiles to be roughly 
constant with radius and slightly tangentially biased on av- 
erage. We do not find significant differences between the 
profiles for cored and NFW models (the reason for this will 
become clear in ti4.3p . Fornax' s anisotropy /3 ^ —0.2 ±0.2, 
while Sculptor and Carina have on average (3 ^ —0.5 ± 0.3. 
For Sextans the anisotropy cannot be determined reliably, 
/3^— 0.3±0.5. These values are compatible with those of 
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Fig. 3. Pdf for the two free parameters characterizing the dark; halo profiles for each dSph galaxy obtained using Schwarzschild 
modeling. The top row shows the pdfs with NFW/Einasto models, the bottom panel those for all cored models explored. The 
contours show the 1 and 2cr confidence levels (the 3a contour is not shown to avoid crowding the image). 



I Walker et al.l ([2003), which were derived using the spheri- 
cal Jeans equation assuming a constant velocity anisotropy 
profile. 

4.2. Bayesian comparison of the models 

We compute the evidence relative to the NFW using Eq. ^ 
by integrating over the parameters (in our case the scale 
radius and the mass) the pdfs shown in Fig. [3] We do 
this for each dwarf galaxy and for all the models listed in 
Table 2) The different Bayes factors are shown in Fig. [S] 
Each set of bars shows the Bayes factors for the given dSph 
galaxy {Bi^NFW,k), while the last set shows the combined 
result {Bi^NFW,comh = Y{k^i,NFW,k)- Wc notc that an 
odds ratio be tween 1:2 till 1:3 is considered "Barely worth 
mentioning" ([JefFrevJ 1 19981 ). and only odds ratios above 
1:10 are considered "strong" evidence. 

For each galaxy there is hardly any evidence for or 
against an Einasto profile (with a' = 0.2, 0.4) compared 
to NFW. This is not unexpected since these profiles are 
quite similar over a large region (|Vera-Ciro et al.l |2013[ ) . 
Also in the case of the combined evidence the NFW and 
Einasto are hard to distinguish. Comparing the NFW or 
Einasto profiles for individual galaxies to the cored mod- 
els, one cannot strongly rule out a particular model. For 
Fornax, Sculptor and Carina, the 7 = 2 models (where the 
transition speed is fast) appear to be less likely, but this 
is not the case for Sextans. However, when we look at the 
combined evidence, i.e. we explore whether all dwarfs are 
embedded in the same halos, such 7 = 2 models are clearly 
disfavored. 

The results for Sc ulptor may be compared to those of 

iBred dels et al] ()2012D . In that paper, the authors found 
that the maximum likelihood value for the central slope of 
the density profile corresponded to a cored model. Since the 



evidence is the integral of the pdf, and not directly related 
to the maximum likelihood (except for a Gaussian distribu- 
tion), we should not be surprised to find a slightly stronger 
evidence for the NFW case here. In any case, the differ- 
ences between the models are minor as shown graphically 
in Fig. m and the evidence and the maximum likelihood 
(marginally) favoring different models can be attributed 
simply to not being able to distinguish amongst these. 



4.3. A robust slope measurement 

We now inspect in more detail the shape of the mass dis- 
tributions found for the various best fitting models. We 
are interested in exploring why the differences between the 
various models as small as apparent in Fig. |4l 

The top row of Fig. [5] shows the enclosed dark matter 
mass for the best fit models (indicated by the solid dots in 
Fig. [3]) for each galaxy separately. We use the same color 
coding as in Fig. [31 and also include the stellar mass in 
black. The red-dashed vertical lines denote r_3, the radius 
at which the light density profile has a logarithmic slope 
of —3, while the black line indicates ri/2- This remarkable 
figure shows that for each galaxy there is a region where the 
mass distributions are truly almost indistinguishable from 
one another. The different profiles, each characterized by 
its own functional form, scale radius Ts and mass Mikpc, 
conspire to produce a unique mass distribution. This re- 
gion extends from slightly below r_3 to approximately the 
location of the outermost data point (see bottom panel). 
Here M{r) ex r^, where x ranges from 1.65 for Fnx, to 1.9 
for Sextans. 

In the middle row of Fig. [HI we plot the logarithmic slope 
of the dark halo density distribution, where the black line 
denotes the stellar density. Near the position where the 
logslope of the stellar density is —3, all the best fit dark 
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Fig. 4. Similar to Fig. 2, except now we show the different best fit models for the various density profiles explored, which are 
indicated with different colors (the color scheme is the same as in Fig. [S} . 



matter density profiles seem to reach a similar logslope, al- 
though the value of the slope varies from galaxy to galaxy. 
The radius where the logslopes coincide lies, as expected, 
inside the region where the mass distribution is well de- 
termined, since both quantities are related through deriva- 
tives. 

To illustrate the distribution of the kinematic sample 
with respect to the light, we plot in the bottom row of Fig. 
[5] the cumulative 2d radial distribution of the kinematic 
data in black. The cumulative 2d radial distribution for 
the light is plotted as the red histogram. All kinematic 
datasets are more concentrated than the light, but no clear 
trend is visible between the distribution of the kinematic 



sample with respect to the light, and the exact location 
where the logslope of most accurately determined. 

The existence of a finite region where the mass is more 
accurately determined has also been observed in the lit- 
erature in works using MCMC in combination with Jeans 
modeling. F 'or example , it is v isible in e.g . the right pane l 
of Fig 1. in 'Wolf et al.' ('2010'), Fig. 18 in lWalkeil (|2012b[ ). 
and Fig. 10 in Jardel et al. ( 2013 ) for Draco, in the case of 
a non-parametric density distribution with Schwarzschild 
models. 

The analysis of lWolf et al.l ()2010f ) used the light weighted 
average of the velocity dispersion to relate the radius at 
which the logslope of the light is —3, or the half light ra- 
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Fig. 5. Evidences for all models listed in Table [J] relative to the NFW case. The last column shows the combined evidence for 
all galaxies together, and shows that the core23 and core24 are strongly disfavored. 
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Fig. 6. Top row: Enclosed mass as a function of radius for the different dark matter density profiles, with the stellar component 
in black. Middle row: Logarithmic density slope as a function or radius, where the black curve corresponds again to the stellar 
component. The red dashed line indicates r_3, the radius at which the light profile has a logarithmic slope of —3, while the black 
line indicates ri/2, the radius at which half of the stellar mass in enclosed (in 3d). Bottom row: Cumulative density distribution 
of the (2d) radial distribution of the data (black), and the light (red) showing the kinematic data is sampled more concentrated 
towards the center. 



dius, to the point where the mass is accurately (being in- 
dependent on the anisotropy) and precisely (showing the 
least uncertainty) determined. Our findings go beyond this 
result. They suggest that whatever dynamical model or 
method is explored, there is a better set of parameters to 



describe the mass distribution of dSph galaxies. Let r_3 
be the radius at which the logslope of the (3d) light distri- 
bution is —3. Since the mass is accurately determined in 
this region, a natural parameter would be M_3 = M(r_3). 
And since also the logslope at this radius is accurately de- 
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Fig. 7. Similar to Fig. [S] except now using M_3 and r_3 as parameters. Note that the contours for the NFW cannot go beyond 
K(r_3) > -1. 



termined, the next parameter should be k_3 = g^^^|r=r_3- 
For any general model, if the values of /3 and 7 are fixed, 
this effectively makes r^ a function of k_3. 

Fig.[7]shows the pdf for the M_3 and k_3 parameters for 
both the NFW and corelS models for each galaxy, assum- 
ing a flat prior on these parameters in the domain shown 
in this figure, except for the NFW profile which we limit 
to K_3 = —1.05, since for k_3 > — 1 the scale radius is 
unphysical. As can be seen from the pdf while there is still 
uncertainty associated to the logslope at this radius, the 
value is nearly model independent and therefore we believe 
this value to be accurate, especially for Fornax and Sculp- 
tor (k_3 = —1.4 ±0.15 and k_3 = —1.3 ±0.12 respectively 
for the NFW model). Note also that some uncertainties 
might arise because the kinematics are not sampled exactly 
according to the light. 

These results also help us understand why we found that 
the scale radii of the best fitting NFW profiles always to be 
larger than those of the cored models (see Fig. [3]) . For the 
NFW, we have 



K[r) 



dlogp(r) 



-2r 



1, 



d log r r + Ts 

which can be easily solved for Ts 

K + 3 



^s.nfw 



K+ 1 



(8) 



(9) 



A similar solution can be found for the other parametric 
models, for instance the 7/? model gives: 



'■s,7/3 



P + K 



-1/7 



(10) 



If we now require that the slopes are the same at r_3 for 
the NFW and corelB models, we find 



r-s.coreia K + 1 ' 



(11) 



which is > 1 for k < — 1, explaining why the cored profiles 
have smaller scale radii than the NFW profile, i.e. to get 
the same logslope at the same location for the cored models, 
their scale radius needs to be smaller than that of the NFW 
profile. A similar result holds for the other cored models. 



5. Conclusions 

In this paper we have presented a comparison of dynam- 
ical models using different dark matter profiles for four 
dwarf spheroidal galaxies in the Local Group, namely For- 
nax, Sculptor, Carina and Sextans. The model comparison 
was done using Bayesian evidence. We have found that no 
particular model is significantly preferred, and that all four 
dwarf spheroidals are compatible with either NFW/Einasto 
or any of the explored cored profiles. Only Sextans shows a 
slight preference for cored models, but not with high odds. 
Nonetheless, we find that it is very unlikely that all four 
dwarf spheroidals are each embedded in a cored dark mat- 
ter halo of the form pom oc 1/(1 + r^)''/^, with /3 = 3, 4. 

Our best fit models however, conspire to produce the 
same mass distribution over a relatively large range in radii, 
from r_3 up to the last measured data point (which is of- 
ten close to the nominal tidal radius obtained from fitting 
the light profile). This M{r) ^ r^, with x = 1.65- 1.9, 
is similar to that suggested by I Walker et al.l ()2009br ) al- 
beit with a slightly steeper exponent. Another (related) 
quantity that is robustly determined and independent of 
the assumed dark matter density profile, is the logslope of 
the density distribution at r_3. We find for the dwarfs 
in our sample, that this slope ranges from k_3 ^ —1.4 at 
r_3 = 0.96 kpc for Fornax, to k_3 ^ —1.1 at r_3 = 0.98 
kpc for Sextans. 

These findi n gs can be seen as an extension of the results 
of IWolf et al.l ( 2010r ). who showed that the mass at r_3 
can be determined very accurately in a model independent 
fashion. These authors demonstrated that this result might 
be understood from the Jeans equation. Although we do 
not have yet a solid mathematical explanation for our new 
findings, we suspect that this might be obtained using the 
virial theorem, which is effectively another, yet independent 
moment of the collisionless Boltzmann equation. 

In the near future, we will apply Schwarzschild modeling 
to the same data but instead of the moments, we will use the 
discrete individual measurements directly. This approach 
should allows us to get the most out of the data, since no 
information is lost. When binning, one loses spatial reso- 
lution, but also the higher moments of the line-of-sight ve- 
locity distribution are not included in the fitting procedure 
because of their large and asymmetric errors. Furthermore, 
the use of the full line-of-sight velocity distribution should 
improve the precision of the anisotropy profile, which may 
be an interesting quantity to discriminate formation sce- 
narios. 
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This moments-to-discrete modeling step must be carried 
out before deciding if and how much more data is needed 
to discriminate among various dark matter density profiles. 
Nonetheless, we have learned here that the functional form 
of the mass distribution may be determined over a large 
distance range, even when only a few hundred velocity mea- 
surements are available (as in the case of Sextans). How- 
ever, the uncertainty on the value of the exact slope of the 
density profile at e.g. r_3 is driven by the sample size. 

An obvious next step is to establish if the subhalos ex- 
tracted from cosmological simulations have the right char- 
acteristics to host the dSph of the Milky Way, now that 
not only the mass, but also its functional form (1st and 
2nd derivatives) , of their dark halos have been determined 
reliably. 
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